Mid-air collisions enhance saltation 
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Here we address the old question in Aeolian particle transport about the importance of mid-air 
collisions. We find that surprisingly these collisions do enhance the overall flux substantially. The 
effect depends strongly on restitution coefiicient and wind speed. We can explain this observation 
as a consequence of a soft bed of grains which floats above the surface and reflects the highest flying 
particles. We make the unexpected observation that the flux is maximized for an intermediate 
restitution coefiicient of about 0.65, which is comparable to the values experimentally measured for 
collisions between sand grains. 



PACS numbers: 45.70.Mg,47.27.-i,47.55.Kf, 83.80. Hj 

Would a sandstorm be stronger if the sand grains in 
the air would not collide against each other? This ques- 
tion has puzzled practitioners and theoreticians alike in 
the past. Models for Aeolian sand flux in Refs. THB] be- 
come certainly much simpler if such mid-air collisions are 
neglected, but does this approximation underestimate or 
overestimate the value of the calculated saturated flux? 
As opposed to experiments, the direct computer simu- 
lation of saltation, the main Aeolian transport process, 
offers the possibility of switching on or off the collisions 
between particles or of modifying at will the collision pa- 
rameters, like for instance the coefficient of restitution, 
allowing in this way to determine for the first time with 
precision the role of mid-air collisions during the saltation 
process. 

We discover that mid-air collisions are in fact the key 
ingredient to understand the relation between different 
concepts that have been put forward in the past like the 
splash [71 IH] , the soft bed [SJ [TU] , and the distinction be- 
tween saltons and reptons [TTJ [T^] . In saltation, particles 
are ejected from the granular bed in a splash, produced 
by the impact of fast particles called saltons (yellow tra- 
jectory in Fig. [T]). These saltons must have the necessary 
kinetic energy to assure that, despite the substantial dis- 
sipation in the ground, some ejected particles can again 
fly high enough. After all, only high-flyers can acquire 
sufficient acceleration to become again saltons, because 
the wind velocity at the ground is zero, increasing loga- 
rithmically with height. Our detailed study reveals the 
following picture: Due to the irregularities of the surface 
a splash produces in fact three types of moving particles 
(see Fig. [I]) : Many (green) creepers that do not leave the 
bed, many (red) leapers making small jumps remaining 
in regions of small wind velocities and that, therefore, 
cannot produce a new splash, and very few saltons (yel- 
low) which fly higher up. Only the saltons assure that the 
saltation process is sustained. Both creepers and leapers 
(reptons) contribute considerably to the global sand flux 




FIG. 1. (color online). Typical splash obtained numeri- 
cally with our model in 3D. The impinging particle (yellow) 
bounces after ejecting other particles from the bed (red and 
green). While the green particles essentially only move on the 
ground, the red ones are lifted and dragged by the wind. 



but play a very different role in what follows. 

Can one really make such a sharp distinction between 
leapers and saltons? After all, the collision process is ran- 
dom and should yield a continuous distribution of ejec- 
tion velocities. The experimental observations say yes: 
Already Mitha et al. [13] observed that the splash distri- 
bution is bimodal as later reproduced numerically by An- 
derson and Haff [14], i.e., that there is a broad spectrum 
of slow particles and a small peak of fast ones, and the 
analysis of these experimental observations led Andreotti 
[TTj to coin the terms saltons and reptons. What makes 
these two types of particles different are the mid- air colli- 
sions. In Fig. [2] we see the trajectory of a typical salton 
simulated with the method described below in in three 
dimension: It makes several jumps without touching the 
ground, rebouncing upwards each time due to a collision 
with a (red) leaper. The probability for such a collision is 
reasonably high because the leapers form a rather dense 
layer which is precisely the soft bed described earlier in 
Ref. [TD]. Consequently, the salton stays longer in areas 
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FIG. 2. (color online). Simulated trajectory of a salton in 
three dimensional saltation. The (yellow, upper trajectory) 
salton is kept in the air by colliding against (red, lower trajec- 
tories) particles from the soft bed with e — 0.7 and for Shields 
number 8 — 0.9. For clarity only the relevant particle trajec- 
tories are shown and the ground is schematically represented 
by a flat plane, (details in the text). 



of strong wind and less time close to the ground, where 
the wind drag acceleration is weaker. This explains why 
the saltons acquire so much energy and can sustain the 
saltation process. Summarizing, mid-air collisions in the 
soft bed are the crucial mechanism that differentiates the 
saltons, essentially doubling the saturated sand flux as we 
will show in the following. 

The Discrete Elements Method (DEM) based on 
particle-particle and particle-wind interactions allowing 
to study quantitatively, the Aeolian transport of granu- 
lar systems [15] . It has been used to calculate the onset of 
saltation and confirm the existence of a jump in the total 
flux in two dimensions [16 . Here, we apply a three di- 
mensional DEM scheme to investigate the role of mid-air 
collisions in saltation. 

Every sand grain is represented by a hard sphere of 
average diameter Dmean- The collisions between parti- 
cles are described in detail in the Supplemental Material 
|17j . We distinguish collisions between particles in the 
bed with restitution coefficient ebed — 0.7 and between 
flying particles (mid-air collisions) with restitution coef- 
ficient e. Gravity acts in vertical direction (y-direction) 
and a logarithmic wind velocity profile u{y) is imposed 
in the horizontal direction (x-direction) , 
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where yo = Dmean/'iO is the roughness of the bed, the 
bed height, k — 0.4 the von Karman constant, and u■^, 
the wind shear velocity. The Shields number defined as 
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FIG. 3. (color online). Concentration and saturated flux, (a) 
Concentration profile of particles as a function of the height 
y for 6 — 0.9 in the absence (red, lower curve) and the pres- 
ence (yellow, upper curve) of mid-air collisions (e = 0.7). (6) 
The relation between the saturated flux and the restitution 
coefficient, for three different Shields numbers exhibits a peak 
around e — 0.65 for higher shear velocities. The arrows show 
the flux without mid-air collisions. 



the grains is also explained in the Supplemental Material 
We consider particles of average diameter 

Dmean = 2 X lO"'* m, sizC dispcrsiou ao = O.WDmean, 

and density ps — 2650 kg/m"^, in a three-dimensional 
wind channel of dimension (700 x 75 x 7.5)D^gQ„ with 
a reflective upper boundary, placed suficiently high to 
avoid any particle colliding against it, and periodic 
boundaries in the other directions. For the fluid density 
we considered ps = 1.174 kg/m'^. In fact, no collision 
with the upper boundary ever occured in our simulations. 
The lower boundary at y = 0, representing the deep 
ground, is strongly dissipative with a fixed restitution 
coefhcient of e^, = 0.5. We consider a bed of 12 particle 
layers to suppress the reflection of shock waves on the 
lower boundary due to the finite depth [TMSU] . In the 
beginning of the simulation, a few particles are dropped 
at random positions to trigger saltation. 

The dimensionless flux in the direction of the wind is 
defined as. 
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is the pertinent dimensionless parameter that controls 
the wind velocity, where s = ps/ Pw the ratio between 
the grain and fluid density, and g the norm of the gravita- 
tional acceleration. The feedback procedure that extracts 
the momentum from the wind due to the acceleration of 



where A — (75 x 7.5)I?^eQ„ is the area of the bottom 
of the channel, vf and rrii are, respectively, the velocity 
along the horizontal direction and the mass of the par- 
ticle z, and D = psy/{s — ^)gD'^^^^ is a normalization 
constant. The saturated flux, which is the average flux 
in the stationary state, does not change with the num- 
ber of particles in the system. The granular temperature 
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FIG. 4. (color online) The temperature (a) and the flux (6) 
profiles for e = 0.7, 1.0, and without mid-air collisions, for 
9 — 0.9. At every height in (a), the granular temperature for 
e = 1.0 is larger than e = 0.7. The flux profiles in (b) confirm 
the higher fiux for e = 0.7. The fiux profile is defined at each 
height as the product of the concentration and the average 
velocity in each slice, (c) The dependence of the maximum 
temperature on the Shields number indicates a linear and a 
sublinear growth for the conservative and dissipative case, 
respectively. 



defined as 



N{y) 



N{y) 
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quantifies the fluctuations around the mean velocity 

The wind channel is divided along the y-direction in 
horizontal shces of volume (2.5 x 75 x 1)D^^^^ to calcu- 
late the profiles of particle concentration, average particle 
velocity, alignment, and granular temperature. The par- 
ticle concentration is the ratio between the volume of par- 
ticles and the total volume. The flux profile is obtained 
for each slice from the product of the concentration and 
the average particle velocity. 

Figure |3]d shows the dependence of the saturated di- 
mensionless flux (Eq. (3)) on e for different Shields num- 
bers. For > 0.44, the saturated flux is higher for 
e — 0.65 than for e = 1.0. The maximum at e = 0.65 
increases substantially with the Shields number, i.e., the 
stronger the wind, the higher the peak. This is also con- 
firmed by the flux profile in Fig. 4b, where the area 
below the curve, which corresponds to the total flux, 
is larger under dissipation. Interestingly, this optimal 
restitution coefficient is comparable to the values experi- 
mentally measured for collisions between grains of quartz 
sand 21 . Below we will explain this maximum in the 
flow at a specific value of the restitution coefficient as 
the result of the competition between two different mech- 
anisms: the decreasing alignment of trajectories and the 
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FIG. 5. (color online) (a) The angle variance between par- 
ticle trajectories for e = 0.7, e — 1.0, and without mid-air 
collisions, for 6 — 0.9. The inset shows the alignment a for 
the three cases. They confirm the weaker alignment in the 
simulations with elastic collisions. (6) Also for = 0.9, the 
probability distribution for the maximum height, which can 
be fitted by the superposition of two exponential decays (solid 
(red) line): one for the leapers (dash-dotted (green) line) 
[0.3ea;p(— 0.12a:;)] and another for the saltons (dash-double- 
dotted (purple) line) [0.0094ea:p(-0.0182a;)], truncated for 
y > 220 possibly due to finite-size effects. For the sake of 
comparison, we include results without mid-air collisions (blue 
triangles), fitted by a single exponential decay (dashed (yel- 
low) line) [0.072e2;p(— 0.05172;)]. In the inset, we show the 
relation between the fiying time and the maximum height, 
where the horizontal dashed line delimits the soft bed. 



increasing uplift of particles with e. 

In the absence of mid-air collisions, all particles are 
reptons (either creepers or leapers) and they follow 
stretched parabolic trajectories as typically assumed in 
saltation models [1]. With mid-air collisions, saltons 
emerge and they are sustained at higher positions, where 
the wind is stronger, by multiple collisions with leapers 
(as shown in Fig. 2). The number of collisions, and con- 
sequently the flying time of a salton, strongly depends on 
the concentration of leapers. Figure|3^ compares the con- 
centration profiles for 9 = 0.90 with and without mid-air 
collisions. The red (dark) curve obtained without colli- 
sions decays similarly as the one proposed by Jenkins and 
co-workers [32 [23] . In the presence of mid-air collisions 
(yellow (light) curve), particles are uplifted and, above 
a certain height, the concentration of particles is much 
larger than in the case without collisions. 

Simultaneously, another important mechanism is ac- 
tivated. While, without dissipation, particles rebound 
randomly in the air, dissipation tends colliding particles 
to align their trajectories [24|- In the Aeolian transport, 
this symmetry breaking is expected to occur in the direc- 
tion of the wind, affecting the overall flux. To investigate 
this alignment, we measure, for the dissipative (e = 0.7) 
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FIG. 6. (color online) Contribution to the total saturated flux 
of sheet flow (lower curve) , the leapers in the soft bed (middle 
curve), and the saltons above the soft bed (upper curve) for 
e — 0.7 as function of the Shields number. The green area 
(lower region) shows the sheet flow in the particle bed y < ho. 



and the conservative (e = 1.0) cases, the granular tem- 
perature T(jj), the inverse of the temperature anisotropy 
given by a{y) =< w^(y) > /T{y), and the variance of 
the velocity angle with respect to the wind direction, 
i.e., (T^ =< > — < a >^, with a — a,Tctaii{vy/vx). 
The granular temperature profile in Fig. [4^ displays for 
e ~ 1.0 at every height a larger temperature for e — 1.0 
as compared to e = 0.7 with a peak at a certain height. 
As a reference, at zero temperature, all particles would 
flow aligned with the same velocity. Figure [3]: shows the 
dependence of the maximum temperature on the Shields 
number 9 for both cases. The maximum temperature 
grows with 9^^'^ (linear in it*) for e = 0.7, and increases 
linearly with 9 (quadratically with u*) for e = 1.0, be- 
ing always larger in the conservative case. This higher 
temperature could either be due to a larger dispersion in 
the magnitude or in the direction of the particle veloc- 
ities. To distinguish these two possibilities we plot, in 
Fig. [5^, the variance of the velocity angle and the tem- 
perature anisotropy (main plot and inset, respectively). 
These data quantitatively confirm that the lower temper- 
ature in the inelastic case is mainly due to the alignment 
of particle trajectories. Clearly, the alignment of the par- 
ticle trajectories, which tends to enhance the overall flux, 
decreases with e. 

In the absence of mid-air collisions, the distribution 
of maximum heights can be approximated by a Poisson 
process and, therefore, well described by an exponential 
decay, as shown in Fig. [5]d (dashed (black) line). With 
mid-air collisions, however, the distribution can only be 
described by the superposition of two exponentials. The 
first exponential dependence (dash-dotted (green) line) 
corresponds to the leapers, while the second one cor- 
responds to the saltons, which typically fly above the 



leapers (dash-dotted line). We can now quantitatively 
define the upper bound of the soft bed as the height above 
which the saltons are in majority. From the trajectories 
of individual particles we can also relate the flying time 
to the maximum height (see inset in Fig. ^jp) ; we observe 
that saltons stay much longer in the air than leapers. 

Three major types of transport contribute to the total 
flux in saltation: the sheet flow, the transport of leapers 
(in the soft bed), and of saltons (above the soft bed). 
Once identified the limits of the soft bed, for each Shields 
number, we can compute the contribution of each mech- 
anism to the total flux, as shown in Fig. [6j The sheet 
flow, computed from the mass transport in the particle 
bed (y < /iq), results from the wind shear stress and the 
creep of particles on the surface. The contribution of the 
leapers and saltons is obtained from the flux in and above 
the soft bed, respectively. As observed in the figure, the 
relative contribution of the saltons and leapers signifi- 
cantly increases with the Shields number in comparison 
to the sheet flow. Figure |6] also shows that the saltons 
contribute the most to the total flux. 

In the saturated flux, we also reproduced the discon- 
tinuity at the onset of saltation reported for 2D in Ref. 
[16j . We also studied the impact of mid-air collisions 
in two dimensions. The decrease in the spatial dimen- 
sion enhances the relevance of mid-air collisions and the 
consequences discussed here are even more pronounced 
However, as expected, the curves for the saturated flux 
in 2D and 3D overlap if the Shields number in 2D is 
rescaled by an appropriate factor that takes into account 
the width of the system. 

Summarizing, the contribution of mid-air collisions 
cannot be neglected in saltation. The absence of colli- 
sions strongly underestimate the total mass flux. The 
saltons contribute significantly to the flux enhancement 
by acquiring large momentum from the wind and using 
it to eject more particles through more intense splashs. 
We provide a new picture of Aeolian saltation where the 
competition between uplift due to mid-air collisions and 
alignment due to inelasticity optimizes the mass trans- 
port in the presence of dissipative collisions. Interest- 
ingly, the position of the maximum corresponds to the 
restitution coefficient typically observed in granular col- 
lisions. 

Our results are crucial for future studies in Aeolian 
saltation, since they provide qualitative and quantita- 
tive information about the influence of mid-air collisions, 
which should be considered in future modeling. Future 
work might include the aerodynamic lift of the particles, 
turbulent wind speed fluctuations, and the electrostatic 
interaction between particles giving a contribution of long 
range forces. 
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PARTICE DYNAMICS 



Trajectories and velocities of the particles are obtained 
itcrativcly by solving the Newton equations of motion, 
through the velocity-Stormer-Verlet scheme [1]. The con- 
tact between particles is described by the spring dashpot 
potential, which has an elastic and a dissipative contri- 
bution acting on each particle. When two particles i and 
j overlap (i.e. when their distance is smaller than the 
sum of their radii) an elastic force is applied, 
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l/2{di 



(1) 



where A: = 0.5 is a spring constant, di and dj are the 
diameters, rrii is the mass of particle i, and rjj points 
from particle i to j. A dissipation force is also applied 
accounting for the inelasticity of the collision. 



^diss ~ I'^iji 



(2) 



where v^j = v^ — Vj is the relative velocity and 7 is the 
dissipation coefficient. The coefficient of restitution e is 
given by the ratio between the absolute velocities after 
and before the collision. For particle-wall collisions, the 
same forces act as if particle i would collide with another 
particle of diameter [2]. For simplicity, friction and 
rotation are neglected. 



II. THE WIND PROFILE 

A logarithmic velocity field mimics the wind profile in 
the absence of saltating grains in the a;-direction with 



y-ho 

u{y) = — In , 



(3) 



where yo = Dmean/30 is the roughness of the bed with 

Dmean the mean diameter of the particles, /iq the bed 
height, K = 0.4 the von Karman constant, and the 
wind shear velocity. In the presence of grains, the wind 
strength is substantially reduced due to the momentum 
transferred to the grains [3] . The grain stress Tg (y) quan- 
tifies the average horizontal force per unit area / that the 
wind applies on the grains above y [4], i.e.. 



(y) = / f{y')dy'. 



(4) 



where / is the average horizontal force per unit volume. 
The modified wind shear velocity Ux(t/) is the fluid stress 
left at the height y after the momentum transfer. 



(5) 



where pw is the air density. To obtain the modified wind 
profile [5], we solve. 



du 
dy 



Urjy) 

ny 



(6) 



considering Ur{y) constant within an interval dy. 

The numerical solution is achieved iteratively. The po- 
sition of the bed surface is used as starting point of the 
integration of the wind profile. However, the particle 
splash and the sheet flow dynamically affect the shape 
of the bed surface. Consequently, we need to find ho for 
every time step. We use a criterion based on the wind 
values to define an approximate position for the particle 
surface. High density areas strongly reduce wind veloci- 
ties. If the calculated velocity Vi in area is below O.lu,, 
we assume that this area contains the particle bed and 
the velocity is set to zero. This means that is cho- 
sen to be the point where the calculated velocity exceeds 
O.lu*. 

The wind drag is the only external force applied to the 
particles along the x-direction. 



Fd = ^PwCdVr'Vr, 



(7) 



where pu, is the air density and Vr — v — u is the velocity 
difference between particle and wind, with u,. = |vr|. 
The drag coefficient proposed by Cheng [6] is suited 
to model grains with irregular, natural, shapes, and is 
given by, 

3/2 

-, 7-, Pw'^r^mean /c\ 

1 ,Re = , (8) 
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where p. = 1.8702 x 10 ^ kg/ (m.s) is the dynamic viscos- 
ity and Re, the Reynolds number. 
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